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We report on the current status of the Golem project which aims at the construction of a general one-loop 
evaluator for matrix elements. We construct the one-loop matrix elements from Feynman diagrams in a highly 
automated way and provide a library for the reduction and numerically stable evaluation of the tensor integrals 
involved in this approach. Furthermore, we present applications to physics processes relevant for the LHC. 
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1. Introduction 

Precise predictions for signals and backgrounds 
at high energy colliders are required to be able to 
identify and study New Physics in current exper- 
iments. Processes with multi-particle final states 
will play a significant role in this task, in par- 
ticular at the LHC, which emphasises the need 
for predictions at next-to-leading order (NLO) for 
such processes. 

One of the challenges of NLO calculations is 
the numerically stable evaluation and integra- 
tion of the virtual corrections. Both the de- 
velopment of new unitarity based methods and 
the improvement of traditional methods based 
on Feynman diagrams have led to important 
new NLO results, in particular to physical pre- 
dictions for 2 — >■ 4 processes at hadron collid- 
ers |l|2|3|4|5|6|7|8|9|10|llj . A recent overview can 



also be found in |12I13| . 

As these calculations are technically involved, 
it is important to spend some effort on two fea- 
tures: automatisation and modularity. The lat- 
ter allows for a combination of automated pro- 
grams to calculate the virtual corrections with 
existing tools for tree level matrix element gen- 
eration, phase space integration and automated 
tools for the subtraction of soft/collinear singu- 
larities from real radiation via a standard inter- 
face [H]. 

The Golem approach to the calculation of one- 
loop amplitudes puts particular emphasis on au- 
tomatisation, as the setup can be used for mass- 
less as well as massive particles and all types of 
spin structures. The Feynman diagrammatic ap- 
proach implies that the rational parts of the am- 
plitudes can be evaluated at no extra cost. The 



reduction of the tensor integrals is done semi- 
numerically, thus avoiding large algebraic expres- 
sions. A special choice of the integral basis pre- 
serves numerical stability in phase space regions 
with small Gram determinants. More details 
about the method are given in sections [2] and |3l 
The Golem method in combination with au- 
tomated tools for the real radiation part has 
been used to obtain results for e.g. qq — > bbbb, 
pp — >■ ZZ + jet and pp — ?> graviton + jet, which 
will be described in sections [5] and [6] 

2. Overview of the Algorithm 

For the consistent description of a cross-section 
for a, 2 —i' N particle process at NLO one needs 
to include the square of the tree level matrix ele- 
ment A^' Aj^j' , the real emission of an additional 
massless particle A]^'^^Aj^j'^-^ and the virtual cor- 
rections consisting of the interference term of one- 
loop and tree-level diagrams of the 2 — > iV pro- 
cess, -'^ /v ^ J\[ I "^-C 

The Golem project focuses on the efficient and 
numerically stable computation of the virtual cor- 
rections, as the other contributions can be added 
by interfacing with existing tools through a well 
defined generic interface [MJ. The one- loop am- 
plitude A]^' is represented as the helicity projec- 
tions of a sum of Feynman diagrams. Each dia- 
gram Qa contributes to one or more colour struc- 
tures and can be written as a sum over tensor 
integrals: 
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In our standard approach we use a Lorentz in- 
variant decomposition of the tensor integrals to 
separate the tensor structure from the associated 
integral form factors. The fully contracted tensor 
structure is expressed in terms of spinor products 
and the constants of the physics model, whereas 
the Feynman parameter integrals are expressed 
in terms of form factors which can either be re- 
duced further or be evaluated numerically. An 



implementation of these form factors is available 
in the library golem95. |15I17) 

As an alternative approach one provides the 
numerator function 



p=0 



■Ql^r, 



(2) 



where q — q — n , as an input to some reduction 
at the integrand level 18 19 20 21 22.:, where Afc 
is evaluated for different numerical values of /z^ 
and the four dimensional, complex vector q. 

3. Reduction of the Tensor Integrals 

The reduction of one-loop integrals follows an 
algorithm which has been developed in [15116) . 
valid for an arbitrary number of legs and both 
massive and massless particles. The ultraviolet 
and infrared divergences are regulated by dimen- 
sional regularisation. 

This algorithm has been implemented 
as a Fortran 95 library, golem95, for 
integrals with up to six external mo- 
menta [17]. The library is publicly available at 
lappweb . in2p3 . f r/lapth/Golem/golem95 . html, 
together with detailed documentation and exam- 
ple programs. 

The algebraic reduction of higher rank four- 
point and three-point functions to expressions 
containing only scalar integrals necessarily leads 
to inverse Gram determinants appearing in the 
coefficients of those scalar integrals. These de- 
terminants can become arbitrarily small in some 
phase space regions and can therefore hamper a 
numerically stable evaluation of the one-loop am- 
plitude upon phase space integration. In order to 
deal with such situations in an automated way, 
our algorithm is able to avoid the occurrence of 
inverse Gram determinants completely. This is 
done by testing the value of the determinant and, 
in the case when it is smaller than a certain value, 
stopping the reduction before such inverse deter- 
minants are generated. In these cases, the tensor 
integrals, corresponding to integrals with Feyn- 
man parameters in the numerator, are evaluated 
by means of numerical integration. The use of 
one-dimensional integral representations guaran- 
tees a fast and stable evaluation. For integrals 



with internal masses, the option to evaluate the 
tensor integrals numerically prior to reduction, 
in regions where the Gram determinant tends to 
zero, is not yet supported, but is under construc- 
tion. For the IR finite scalar box and triangle 
functions with internal masses, we use the ex- 
pressions from the program DneLDop ^23]. The 
inclusion of an option to allow complex masses is 
underway. 

To summarise, the library contains the follow- 
ing features: 

• all tensor coefficients up to rank six six- 
point integrals, for massless as well as mas- 
sive integrals 

• all scalar boxes, triangles, bubbles, tadpoles 

• boxes in n -f 2 dimensions up to rank three 
(n = 4 — 2e) and in n -I- 4 dimensions up to 
rank one, and triangles in n -t- 2 dimensions 
up to rank one (larger ranks of the higher 
dimensional integrals are not needed in the 
reduction). 

• automated detection of small Gram deter- 
minants and fast (one-dimensional) numer- 
ical integration of tensor integrals to avoid 
the occurrence of small inverse determi- 
nants 

• option to calculate the rational parts only. 



4. Automated Construction 
Loop Amplitudes 



of One- 



For the automated construction of the ma- 
trix elements we have developed a package 
(golem-2.0 [21) combining QGRAF [23^ as a di- 
agram generator with the symbolic manipulation 
program Form j26j for the algebraic simplifica- 
tion of the generated amplitude. The expressions 
for the helicity projections of each diagram con- 
tributing to the process are further processed by 
haggles [27j to produce optimised Fortran 95 
code for a fast numerical evaluation. A large part 
of the infrastructure and the user interface is im- 
plemented as a Python program. 

For the reduction of the tensor integrals 
golem-2 . supports two choices. One can express 



the amplitude in terms of integral form factors as 
defined in 15 and implemented in golem95 [17]. 
Alternatively, the program can generate the nu- 
merator of each Feynman diagram in a form suit- 
able as input for Ssunurai [22 . For each of the 
two options the amount of human intervention is 
limited to setting up an appropriate configuration 
file. 

In order to further improve the applicability 
of the program we are currently implementing an 
interface to FeynRules [28] model files. Interoper- 
ability with existing Monte Carlo programs will 
be achieved by providing an interface conform- 
ing to the Binoth-Lcs Houches Accord [Mj and 
some additional functionality to provide the rele- 
vant information about the matrix element to the 
POWHEG BOX [29 . 

The described features offer a great flexibility 
to the potential user; the achieved performance of 
the matrix element code and its modularity allow 
for a wide range of applications. 

5. Results for qq — > bbhh 

We have used the setup described in Section |4] 
for generating a Fortran 95 implementation of 
the one-loop matrix element for the process qq — > 
hhhh [7], which is one of the channels of the full 
hhhh cross-section at the LHC. 

The QCD induced contribution to the 4&- 
signature constitutes an important background 
for Higgs searches in the MSSM, where for cer- 
tain parameter regions the heavy Higgs boson de- 
cays predominantly into 66 pairs 13 30J. Also 
outside the framework of the MSSM interesting 
physics scenarios such as certain hidden valley 
models lead to collider signatures with many 6- 
quarks [5T] . 

We have ensured the correctness of our one- 
loop matrix element by comparison to a second, 
independent setup using FeynArts [32], Form [26] 
and Maple as a tool chain. Integrated results have 
been obtained by three different setups: 

• the virtual matrix element obtained from 
golem-2 . has been used as a stand-alone 
program to reweight unweighted leading or- 
der events generated by WHIZARD [33 51]. 
All tree like parts and the infrared subtrac- 



tion terms as defined in [3S] using a slic- 
ing parameter [3 6) have been implemented 
in WHIZARD. 

• the virtual matrix element from golem- 2 . 
has been linked into MadEvent 37 together 
with the tree-like parts from MadGraph [3S] 
and the dipoles from MadDipole |39l40j . 

• the third setup is the same as the pre- 
vious but using Samurai [22] in place of 
golem95 [T^ for the reduction of the nu- 
merators of the virtual matrix element. 

The results use LHC kinematics {^/s — 14TeV 
centre of mass energy). As we calculate observ- 
ables for well separated jets, we cluster the final 
state particles by applying the fcr-algorithm [41 
into 6-jets, which then have to pass the set of cuts, 



PT{bj) 



> 30 GeV, 
< 2.5, 



(3) 



l\R{b,,bj) > 0.8. 



We use the parton density functions 
CTEQ6.5 [42] with two- loop running of a^. 
The central renormalisation scale is /io = 

\\jYl,j=i\PT{bj)Y while the factorisation scale 
has been fixed at jjlf = 100 GeV. The initial state 
includes q £ {u,d, s,c}; for the whole calculation 
we use the approximation mi, — and rrit — >■ oo. 

During phase space integration numerical sta- 
bility was monitored by checking the magnitude 
of the X-factor for a given phase space point and 
also the quality of the cancellation of the infrared 
poles. In the worst case about 0.5% of the events 
needed evaluation in quadruple precision. By us- 
ing the setup with Samurai, on top of previous 
tests we could also use several reconstruction tests 
as described in [22] for estimating the numerical 
quality for each phase space point. 

Figure [T] shows the distribution for the invari- 
ant mass of the hardest pair of 6-jets, where the 
jets are ordered by their pT- The error bands 
are obtained by a variation of the renormalisa- 
tion scale fiji = xfiQ for x £ [5; 4]. The dashed 
line corresponds to the leading order result at 
Mfl = Mo- As expected, the plot shows a clear 
reduction of the scale dependence. 




Figure 1 . Distribution of the invariant mass of the 
hardest 6-jet pair in qq — ?> bbbb. The error bands 
represent a variation of the renormalisation scale 
fiR = x^iQ for X £ [i;4]. 



In Figure [2] we examine the precision of 10^ 
phase space points by three different criteria. 
Since the cancellation of the single and double 
poles on their own are not a reliable indicator for 
numerical problems we also consider the local K- 
factor as a heuristic criterion (see [43]). In the 
chosen setup with MadEvent as an integrator and 
Scunurai for the reduction of the amplitude the 
number of unstable points in double precision is 
below one per mil. 

6. Results for pp -^ T^F+jet and pp — > 6"+ jet 

Golem tensor reduction methods have also been 
used to compute fully-differential NLO QCD cor- 
rections to ZZ+jet [44] and Kaluza-Klcin gravi- 
ton-l-jet [45] production at the LHC and Teva- 
tron. At these colliders, graviton-f jet production 
is an important channel for graviton searches and 
ZZ+jet production is an important background 
for Higgs particle and new physics searches. In 
both calculations, we have identified the renor- 
malisation and factorisation scales: fj,R = fip = 
/i. In Figs. |3] and |4] we show the scale depen- 
dence of LHC cross sections for ZZ+jet and 
graviton-|-jet production, respectively. The cen- 
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Figure 2. Criteria used for the determination 
of numerical stability of the virtual matrix ele- 
ment: cancellation of the infrared poles (lower 
cc-axis) and pointwise iiT-factor (upper a:-axis) for 
10^ events. The virtual part has been evaluated 
using Samurai in double precision. 
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tral choice for /i is the Z mass and the trans- 
verse momentum of the graviton p^ , respectively. 
The graviton search selection requires p™"**^ > 500 
GeV and jets are required to satisfy p^ > 50 GeV 
and |?7j| < 4.5. For ZZ+jet production, we re- 
quire pt > 50 GeV for the hardest jet. Following 
Ref. [HI, in Fig. |3]we also show the NLO cross 
section when 2-jet events with a second hardest 
jet with Pt > 50 GeV are vetoed. The jet veto 
clearly reduces the residual scale variation. We 
employed the CS dipole subtraction 35j imple- 
mentations in MadDipole [39 and SHERPA [47] to 
calculate numerical results for the finite real cor- 
rections contribution. For both calculations in- 
ternal cross checks and comparisons have been 
carried out (see |44|45j ). Our ZZ+]et calculation 
has also been validated through a tuned compari- 
son with the independent calculation of Ref. [35] • 
For a fixed phase-space point, the virtual correc- 
tions, obtained using different calculational tech- 
niques, agree at the level of 10^* or better. The 
comparison of full NLO cross sections for the LHC 
and Tevatron, shows agreement at the per mil 
level. Details are reported in Ref. [T^. 



Figure 3. Scale dependence (/x_r, — fip — fJ-) of 
the ZZ+jet cross section at the LHG (-^/s = 14 
TeV) with prjct > 50 GeV for the hardest jet 
in LO (dotted) and NLO (solid). The exclusive 
NLO cross section when a PTjet > 50 GeV veto 
for additional jets is applied is also shown (dot- 
dashed). 



7. Conclusion 

We have given an overview of the recent devel- 
opments in the Golem collaboration whose main 
goals are the automatisation of the computation 
of one-loop matrix elements and the provision of 
precise predictions for LHC observables. As an 
integral part of this work we developed a library 
(golem95) of form factors for the reduction of one- 
loop tensor integrals, which in its latest version 
contains all integrals for one to six propagators 
up to full rank both for massless and massive 
particles in dimensional regularisation. We have 
implemented the Golem method for computing 
one-loop amplitudes in a highly automated frame- 
work (golem-2 . 0) which has been applied to the 
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Figure 4. Scale dependence {fifj = fip = /^) of the 
graviton+jet cross section at the LHC (-y/s = 14 
TeV) in LO (dotted) and NLO (soUd) in a Kaluza- 
Klein model with (5 = 4 extra dimensions and a 
fundamental scale Ms = 5 TeV. Selection cuts 
are described in the main text. 



computation of the quark-induced channel of the 
process qq -^ bbbb. A similar method using the 
Golem tensor reduction has been used to com- 
pute the production of Z-pairs associated with a 
jet and the production of a graviton plus one jet 
at the LHC. 
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